Topological solitons in three-band superconductors with broken time reversal symmetry 
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We show that three-band superconductors with broken time reversal symmetry allow magnetic flux-carrying 
stable topological solitons. They can be induced by fluctuations or quenching the system through a phase 
transition. It can provide an experimental signature of the time reversal symmetry breakdown. 
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Experiments on iron pnictide superconductors suggest the 
existence of more than two relevant superconducting bands 
[1, 2]. The new physics which can appear in these circum- 
stances is the possible superconducting states with sponta- 
neously broken time reversal symmetry (BTRS) as a conse- 
quence of frustration of competing interband Josephson cou- 
plings [2] (other scenario for BTRS state was discussed in 
[3]). BTRS states also attracted much interest earlier in the 
context of unconventional spin-triplet superconducting mod- 
els. There they have a different origin and are described by 
two-component Ginzburg-Landau models [4]. In those cases 
the theory predicts domain walls which pin vortices [4] . It was 
suggested that this can result in formation of experimentally 
observable vortex sheets if (i) a domain wall itself is pinned 
by sample inhomogeneities, or (ii) if a domain is dynamically 
formed inside a current-driven vortex lattice [4]. 

Here we show that a BTRS state in a three-band supercon- 
ductor allows formation of metastable topological solitons. 
Although it is not by any means required to be near Tc for 
these solitons to exist, we use a static three-band Ginzburg- 
Landau (GL) free energy density model : 
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Here, D = V ^ ieA, and ifji = l^/^ile*"^^ are complex 
fields representing the superconducting components. We 
choose to work here with a minimal effective potential V = 
X]i=i 2 3 ^*IV^*P + ^ftl^il^- Although there could be vari- 
ous other terms allowed by symmetry in (1) they are not qual- 
itatively important for the discussion below. For r]ij > 0, 
the Josephson interaction term is minimal for zero phase dif- 
ference, while r]ij < it is minimal for (pi — (pj = tt. 
When the signs of r]ij coefficients are all positive, [we de- 
note it as (+ + +)] the ground state has (fi = (p2 = (fs. 

Similarly in case (H ) one has phase locking pattern 

^1 = V^2 = ^3 + TT. However in cases (+ H — ) and ( ) 

there is a frustration between the phase locking tendencies 
[i.e. one cannot simultaneously satisfy cos{(pi — (pj) = ±1]. 
For example, consider the case a^ = — 1, Pi = 1 and 
r]ij = —1. Without loss of generality lets set cpi = then 
two ground states are possible (p2 = 27r/3, cps = — 27r/3 or 
(p2 = — 27r/3, (ps = 27r/3. Thus in these frustrated cases 
there is Z2 broken symmetry in the system associated with 



complex conjugation of the all tjj fields. The broken Z2 sym- 
metry implies existence of domain walls solutions, which are 
schematically shown on Fig. 1. Note that the frustrated phase 
differences can assume values different from 27rn/3 in case of 
differing effective potentials or Josephson coupling strengths. 
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Figure 1. (Color online) - Schematic representation of various Z2 
domain walls in three -band superconductors with different frustra- 
tions of phase angles, shown by arrows of different colors. Pink line 
schematically shows phase difference between red and green arrow, 
interpolating between the two inequivalent ground states. 



Let us now outline basic properties of the model (1). With- 
out intercomponent Josephson coupling and a^ < 0, its sym- 
metry is [/7(1)]^. Then it allows three kinds of fractional flux 
vortices with logarithmically diverging energy [5] character- 
ized by a phase winding in (i.e. integral over a phase gradi- 
ent around a vortex) Acpi = ^^ \/(pi = 27r. Such a vor- 
tex carries a fraction of magnetic flux quanta (^0), given by 
^- = |^.|2/(|^-^|2^|^2p + |V^3p)<^o. However a bound state 
of three such vortices (i — 1,2,3) has a finite energy. The 
finite-energy bound state is a "composite" vortex which has 
one core singularity where \\\)\\ + \\\)2\ + \^?\ — 0. Around 
this core all three phases have similar winding A(/p^ = 27r. 
Thus it is a logarithmically bound state of fractional vortices 
whose flux adds up to one flux quantum <l>o. In case of non- 
zero Josephson coupling fractional vortices are bound much 
stronger since they interact linearly [5]. 

We show below that the model (1) remarkably has a dif- 
ferent kind of stable topological excitations distinct from vor- 
tices. Note that in two-component superconductors Skyrmion 
and Hopfion topological solitons can be represented as bound 
states of two spatially separated fractional vortices [6]. Like- 
wise we can represent a topological soliton carrying N flux 
quanta (i.e. with each phase winding 2i^N ) in a three compo- 
nent superconductor like a stable bound state of spatially sep- 
arated 3N fractional vortices. Below we will call it ''GL^^^ 
soliton". At first glance, split fractional vortices could not be 
stable in the model (1) because of the strong linear attractive 



interaction between fractional vortices caused by Josephson 
couplings. However we show that such solutions exist as topo- 
logically nontrivial local minima in the energy landscape of 
the model (1). These solutions may also be viewed as combi- 
nations of fractional vortices and closed domain walls. 

Domain walls can form dynamically by a quench, but due to 
its line tension a single Z2 closed domain wall (i.e. a domain 
wall loop) should rapidly collapse. Because of the field gra- 
dients, the superfluid density is suppressed on a domain wall. 
Therefore it can pin vortices. Furthermore at a domain wall 
one has energetically unfavorable values of cosines of phase 
differences cos{(pi — (fj)- Thus Josephson terms immediately 
at the domain wall energetically prefer to split integer flux vor- 
tices into fractional flux vortices since it allows to attain more 
favorable phase difference values in between the split frac- 
tional vortices. (Note that, away from domain walls, Joseph- 
son terms give in contrast attractive interaction between frac- 
tional vortices). We find that if the magnetic field penetration 
length is sufficiently large, then there is a length scale at which 
repulsion between the fractionalized vortices pinned by do- 
main wall counterbalances the domain wall's tension. It thus 
results in a formation of a stable topological soliton made up 
of 3A^ fractional vortices. Thus these topological solitons rep- 
resent a closed Z2 domain wall along which there are N points 
of zeros of each condensate \ijji\. Around each of these zeros 
the phase (pi changes by 27r. The total phase winding around 
the soliton is §Vipidl — ^\/(p2dl = §Vip^dl = 2t:N. 
Therefore it carries N flux quanta. 

Since it is a complicated nonlinear problem, no analytical 
tools are available and thus a conclusive answer if these soli- 
tons are stable could only be obtained numerically. We per- 
formed a numerical study based on energy minimization using 
a Non-Linear Conjugate Gradient algorithm showing the exis- 
tence and stability of the GL^^^ solitons. Technical details of 
numerical calculations are discussed in Appendix . The gen- 
eral tendency which we observed is, that in contrast to most of 
the known topological solitons, they are more stable at higher 
topological charges. In fact we did not find any stable solitons 
for the lowest topological charge corresponding to enclosed 
one quanta of magnetic flux (TV = 1). The lowest topologi- 
cal charge solutions we found carry two flux quanta, and thus 
consist of six fractional vortices residing on a closed domain 
wall. The Fig. 2 shows the TV = 2 soliton in a supercon- 
ductor with two passive bands (thus in this respect, similar to 
the models which are believed to be relevant for iron pnic- 
tide) coupled to an active band. Although it consists of six 
fractional vortices, one of the bands in this example has larger 
density and thus the magnetic field has two pronounced peaks 
near singularities in the main band. This is because the frac- 
tional vortices in that band carry the largest amount of the 
magnetic flux <l>3 = \^2.? /[\i^i? + \^2? + \^2.?]- So the 
magnetic field profile of this soliton resembles a vortex pair. 
We similarly found TV = 2 solitons for superconductor with 
three passive bands and for three active bands which was not 
qualitatively different from the one shown on Fig. 2. 

We find that solutions with larger number of flux quanta 




Figure 2. (Color online) - TV = 2 topological solitons for two 
similar passive bands {ai^^i) — (1,1) with interband coupling 
?7i2 = —3. These bands have Josephson coupling 7713 = 7723 = 1 
to the third band, which is active (a3,/^3) = (—2.5,1). The 
system is type-II with e = 0.07 (we use coupling constant e in 
(1) to parametrize inverse penetration length). The panel A dis- 
plays the magnetic field B. Panels B and C respectively display 
{ipiip2 — '0i'02)/2^ and ('0i'03 — '0i'03)/2i, showing the phase dif- 
ference between two condensates. Second line, shows the densities 
of the different condensates |V^i|^ (D), |V^2|^ (E), |^3|^ (F). The 
third line displays the supercurrent densities associated with each 
condensate | Ji| (G), | J2I (H), | J3I (I). Phase differences on pan- 
els B and C show that there is a closed domain-wall since there 
are two areas with different phase-lockings (blue and red) associated 
with two possible ground states. The solution consists of TV = 2 
vortices which are fractionalized : indeed, the panels D, E and F 
show separated highly asymmetric pairs of singularities of different 
condensates. Note the very complicated geometry of supercurrent 
densities shown on panels G, H and I. 



tend to have ring-like shapes. The Fig. 3 gives an example 
of a solution with TV = 8 flux quanta. Note that this ob- 
ject will have a very distinct magnetic signature which can be 
distinguished by scanning SQUID or Hall or magnetic force 
microscopy. Despite the fact that this object is a bound state 
of 24 fractional vortices, the magnetic field has only 8 pro- 
nounced maxima. They coincide with the position of the 8 
singularities in the band with the largest density. 

The magnetic structure of the soliton always clearly reflects 
the relative densities the bands. When the ground state densi- 
ties in each band are equal, the magnetic field has a uniform 
ring-like geometry as shown on Fig. 4. 

When disparity of the densities in different bands is small 
there is also a family of TV quanta solitons which have 2 TV 
pronounced maxima in the magnetic field. An example with 
TV = 4 is shown on Fig. 5. 

We investigated numerically more than 500 parameter sets 
in three-component BTRS GL models. For all type-II three- 
component BTRS GL models we found stable GL^^^ solitons, 
provided the topological charge was large enough. The solu- 





Figure 3. (Color online) - A^ = 8 quanta soliton for the same param- 
eter set as in Fig. 2 except that e = 0.3 and (as, /^s) = ( — 1-5, 1), 
giving less disparity in the ground state densities (displayed quanti- 
ties are the same as in Fig. 2). The cores of vortices in each bands 
do not coincide. Note the complicated structure of currents in each 
band. 




Figure 4. (Color online) - N — 5 quanta soliton with e = 0.3. With 
three identical passive bands {ai, /3i) = (1,1), with superconductiv- 
ity induced by repulsion r]ij = — 3 between the three condensates. 
Displayed quantities are the same as in Fig. 2. 



tion existed in BTRS states irrespectively of whether bands 
are active or passive and for very different effective potentials 
and interband coupling strengths. It indicates that these soli- 
tons should be rather generic excitations in three-component 
type-II BTRS superconductors. Fig. 6 shows the energy and 
stability of the solitons for different values of the coupling 
constant e (in our parametrization e controls the inverse mag- 
netic field penetration length). It reflects the generic tendency 
which we find, that the solitons are more stable in more type-II 



Figure 5. (Color online) - N = 4 quanta soliton for two similar 
passive bands coupled to a third active band. The parameter set used 
here is the same as in Fig. 3 except (as, /^s) = (—0.5, 1) and e = 
0.2. Displayed quantities are the same as in Fig. 2. 



regimes and also at higher topological charges. 




Figure 6. (Color online) - Energies of the solitons per flux quanta, 
in the units of the energy of a single ordinary vortex (left). When 
the electric charge increases (i.e. the penetration length decreases) 
solitons with smaller N become unstable. The right panel shows 
crossections of the magnetic field for solitons with N G [2,8] 
(double -peak curves). The central curve corresponds to a crossection 
of a regular A^ = 1 vortex. The parameters of the Ginzburg-Landau 
model used here are the same as in Fig. 4, which gives nearly axially- 
symmetric magnetic field. 

Lets us now address the physical observability of these soli- 
tons. First in all the cases which we studied in the model (1), 
the solitons with N flux quanta were more energetically ex- 
pensive than TV isolated one-quanta vortices. However they 
are protected by an energy barrier against decay into ordi- 
nary vortices. Note that because the solitons are obtained 
as solutions of the energy minimization problem, they are 
guaranteed to be stable against infinitesimally small pertur- 
bations. However, since they are more energetic than vortices, 
strong enough perturbation should destabilize them. This sta- 
bility question is addressed numerically in the Appendix . For 
strongly type-II regime the potential barrier can be estimated 
as the energy needed to disconnect the domain wall. For a 
soliton in a three-dimensional sample with phase winding in 



the xy-planc the potential barrier can be estimated as [coher- 
ence length]^ X [sample size in the direction of applied mag- 
netic field] X [condensation energy density] . 

Being more expensive than vortices, these objects cannot 
form as a ground state in low external field [7]. However as 
demonstrated in Fig. 6 they are not much more energetically 
expensive than vortices. In fact the corresponding energy dif- 
ferences can be just a few percent. Thus they can be excited 
by either by (a) thermal fluctuations or (b) by quenching in a 
sample subjected to a magnetic field. To address the scenario 
(b) of possible formation of these solitons in a post-quench 
relaxation, we have to assess "capture basin" of these solu- 
tions (i.e. how large is the area in the free energy landscape 
from which an excited system would relax into the local mini- 



mum corresponding to a soliton. Although studying real post- 
quench relaxation dynamics is beyond the scope of this paper, 
nonetheless we can directly assess the capture basin of the 
solutions from the evolution of the system in our relaxation 
scheme (see also remark [8]). We investigated several hun- 
dreds regimes and found that solitons typically easily form 
when a system is relaxed from various higher energy states. 
This indicates that the capture basin of these solutions is typi- 
cally very large. We find that these defects in fact very easily 
form during a rapid expansion of vortex lattice (which should 
occur when magnetic field is rapidly lowered, or if a system 
is quenched through i^c2)- A typical example is shown on 
Fig. 7. Animations of these processes are available as a sup- 
plementary online material [9]. 




Figure 7. (Color online) - The soliton formation during energy relaxation of an initial state of expanding group of vortices in a circular 
system with open boundary conditions. First line displays the energy density. Second line shows the phase difference between condensates 
(V^r V^2 — V^iV^2)/2^. When domain walls form they separate two inequivalent ground states (blue and red). Third line is the density of the first 
condensate |^i|^. Initial configuration has a high density of 13 vortices in the center. Repulsive type-II interaction makes all vortices move 
away from each other and escape the sample. In the process of energy minimization domain walls and GL^^^ solitons form. Domain wall 
connected to boundaries quickly disappear. The final picture shows the resulting long-living state of a well separated A^ = 4 GL'^^^ soliton 
and a vortex. Parameter set used here is the same as in Fig. 4, with e = 0.4. 



In conclusion, we have shown that BTRS state of a three- 
band superconductor can be detected through its magnetic re- 
sponse. Namely we have demonstrated that in this state the 
system has two kinds of flux carrying topological defects : 
ordinary vortices and also a different kind of topological soli- 
tons. These solitons are only slightly more energetically ex- 
pensive than vortices (in some cases we found the energy dif- 
ference as small as 10~^£^^ where Ey is the energy of a vor- 
tex). They should form during a post-quench relaxation of a 
BTRS superconductor in an external field, since they repre- 
sent local minima with a wide capture basin in the free en- 
ergy landscape. I.e. a system should relax to these local min- 
ima from a wide variety of excited states. Then these solitons 
can be observed in scanning SQUID, Hall, or magnetic force 
microscopy measurements. They can provide an experimen- 
tal signature of possible BTRS states in iron pnictide super- 
conductors. A tendency for vortex pair formation, yielding 



magnetic profile similar to that shown on Fig. 2 was observed 
in Ba(Fei_xCox)2As2, [10] as well as vortex clustering in 
BaFe2-xNixAs2 [11]. These materials have strong pinning 
which can naturally produce disordered vortex states [11], al- 
though a possibility of "type- 1.5" scenario for these vortex 
inhomogeneities was also voiced in [11]. The vortex pairs 
observed in [10] can be discriminated from N = 2 solitons 
(such as that shown on Fig. 2), by quenching the system and 
observing whether or not it forms vortex triangles, squares, 
pentagons etc corresponding to higher- A^ solitons. 
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Appendix A : Finite element energy minimization 

The GL^^^ solitons are local minima of the Ginzburg- 
Landau energy (1). This means that functional minimization 
of (1), from an appropriate initial guess carrying several flux 
quanta, should lead to a GL^^^ soliton (if it exists as a stable 
solution). We consider the two-dimensional problem (1) de- 
fined on the bounded domain (^ C R^, supplemented by a 
'open' boundary conditions on dft. 

Strictly speaking, there is a constraint on dft. This 'open 
constraint' is a particular Neumann boundary condition, such 
that the normal derivative of the fields on the boundary are 
zero. These boundary conditions in fact are a very weak con- 
straint. For this problem one could also apply Robin bound- 
ary conditions on dQ, so that the fields satisfy linear asymp- 
totic behavior (exponential localization). However, we choose 
to apply the 'open' boundary conditions which are less con- 
straining for the problem in question. 'Open' boundary con- 
ditions also imply that topological defects can easily escape 
from the numerical grid, since it would further minimize the 
energy. To prevent this, the numerical grid is chosen to be 
large enough so that the attractive interaction with the bound- 
aries is negligible. The size of the domain is then much larger 
than the typical interaction length scales. Thus in this method 
one has to use large numerical grids, which is computation- 
ally demanding. At the same time the advantage is that it is 
guaranteed that obtained solutions are not boundary pressure 
artifacts. 

The variational problem is defined for numerical compu- 
tation using a finite element formulation provided by the 
Freefem++ library [13]. Discretization within finite element 
formulation is done via a (homogeneous) triangulation over 
ft, based on Delaunay-Voronoi algorithm. Functions are de- 
composed on a continuous piecewise quadratic basis on each 
triangle. The accuracy of such method is controlled through 
the number of triangles, (we typically used 3 ~ 6 x 10^), the 
order of expansion of the basis on each triangle (P2 elements 
being 2nd order polynomial basis on each triangle), and also 
the order of the quadrature formula for the integral on the tri- 
angles. 

Once the problem is mathematically well defined, a numeri- 
cal optimization algorithm is used to solve the variational non- 
linear problem (i.e. to find the minima of J^). We used here 
a Nonlinear Conjugate Gradient method. The algorithm is it- 
erated until relative variation of the norm of the gradient of 
the functional T with respect to all degrees of freedom is less 
than 10-^ 



Initial guess 

As discussed in the paper, N quanta GL^^^ solitons in the 
three-component model are more energetically expensive than 
A^ quanta ordinary vortices. They are local minima of the en- 
ergy functional (1). As a result the initial guess should be 
within the attractive basin of the GL^^^ solitons. Otherwise 



the configuration converges to ordinary vortices which have 
the same total phase winding but cost less energy. We find 
however the attractive basin of the GL^^^ soliton solutions to 
be generally quite large (i.e. the GL*^^^ soliton forms quite eas- 
ily in general). The initial field configuration carrying N flux 
quanta is prepared by using an ansatz which imposes phase 
windings around spatially separated N vortex cores in each 
condensates : 



IV^al =^aY[J^ f 1 + tanh (—(ni(x,y) - ^a) 



A = -— (sin O, — cos 0) 



zG+iAi 



(2) 



where a = 1,2,3 and Ua is the ground state value of each 
superfluid density. The parameter ^a gives the core size while 
6 and 7^ are 



A^^ 



eix,y) = J2^iix,y), 



i=l 



Qi{x,y) = tan" 



1 y-Vi 






N,, 



^{x,y) = ^ni(x,y), 



^i{x,y) = V^(7^^x^)27r(^^^^0^ . 



(3) 



The initial position of a vortex is given by (x^, yi). The func- 
tions Aab = ^b — ^a can be used to initiate a domain wall. 
As an initial guess we generally choose A12 = — A13 = A, 
with A defined as 



{H(r - ro) - 1) 



(4) 



where H(r — tq) is a Heaviside function. Thus in the ini- 
tial guess the domain wall has infinitesimal thickness. It takes 
only a few steps from this initial guess to relax to a true do- 
main wall during the simulations. Consequently, it is entirely 
sufficient to use Heaviside functions for the initial guesses of 
domain walls. Once the initial configuration defined, all de- 
grees of freedom are relaxed simultaneously, within the 'open' 
boundary conditions discussed previously, to obtain highly 
accurate solutions of the Ginzburg-Landau equations. In a 
strongly type-II system when the initial guess was either (a) 
vortices placed on a domain wall or (b) closed domain wall 
surrounding a densely packed group of vortices, the system 
almost always formed GL^^^ solitons. We used also initial 
guesses (c) without any domain walls (A = 0). In that case 
we observed GL^^^ soliton formation, if in the initial states 
vortices were densely packed. This again indicating that the 
GL^^^ solitons in the three component GL model represent 
local minima with wide capture basin in the free energy land- 
scape. 



Figure 7 in the paper shows stages of the energy minimiza- 
tion. Corresponding movies are available as the supplemen- 
tary online material [9] . The main focus of this work is the 
existence of stable static solutions, however the numerical re- 
laxation scheme which we use can give insight into possible 
formation dynamics of these objects. That is, the gauge in- 
variant gradient flow of Ginzburg-Landau free energy can be 
related to the dynamics of Time Dependent Ginzburg-Landau 
equations [8, 14]. Therefore supplementary movies not only 
give information about the size of the capture basin of the lo- 
cal minima associated with the GL^^^ solitons, they also pro- 
vide some insight into possible real dynamics which can lead 
to their formation. 



Appendix B : Stability of the solutions 

The solutions were obtained using an (energy) minimiza- 
tion algorithm, and not by solving the equations of motion. As 
a result, after the convergence (which is carefully controlled), 
the solution is guaranteed to represent (at least) a local min- 
imum of the energy functional (1). Because no symmetry- 
imposing ansatz is used, there are no possible unstable modes 
truncated by symmetry assumptions. Linear stability analy- 
sis consists of applying infinitesimally small perturbation to 
the fields, and investigating the eigenvalue spectrum of the 
(linear) perturbation operator, on the background of a given 
soliton. When the background solution is (meta) stable all in- 
finitesimally small perturbations are positive modes and thus 
can only increase the energy. However a strong perturba- 
tion should cause a decay of a soliton to ordinary vortices 
since these solitons are protected against decay by a finite en- 
ergy barrier. Instead of studying different modes, we double- 
checked the stability numerically by perturbing the solution 
by a random noise. The random noise which is applied to all 
degrees of freedom, is generated as follows 



Re(V^,) = Re(V^,)(o) 
Im(V^,)=Im(V^,)(o) 



PUalJ^a{x^y)i 



i(0) 



Pmax(|A|)/i^(x,y). 



(5) 



Here ^^^ denotes the background solutions, P is a percentage 
giving the relative magnitude of the fluctuation with respect 
to the maximal amplitude of a given field of the background 
solution. /J^{x^y), /j}^{x,y) and ii\{x^y) are (independent) 
random functions of the space G [—1 : 1]. As a result all 
fields initially receive noise whose relative amplitude is P. 
The system is then again relaxed using the same minimiza- 
tion scheme as for constructing the solitons. It is found that 



if the random noise does not exceed a certain threshold, the 
configuration relaxes back to the soliton solution, as can be 
seen from Fig. 8. The noise was gradually increased, finding 
that indeed, sufficiently strong perturbation drives the soliton 
over the barrier, in the energy landscape. Thus leading to its 
decay to ordinary vortex solutions as shown on Fig. 9. The 
precise value of the relative amplitude required to destabilize 
a given soliton, obviously depend on the GL parameters and 
on the number of flux quanta of the solution. 
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Figure 8. (Color online) - Displayed quantities are the same as in Fig. 7 of the paper, namely the energy density, (V^i V^2 — V^iV^2)/2^ and 
I V^i 1^. Initial configuration is a charge 8 GL^^^ soliton shown in Fig. 3. The snapshots show the state of the system a different stages of the 
energy minimization algorithm after the applied perturbation. The initial noise is P = 0.6, which in fact is a very significant perturbation 
where the density fields vary locally up to 60 % of the ground state values, while magnetic field varies up to 60 % of its maximal value. The 
configuration nevertheless relaxes back to the GL^^^ soliton. 




Figure 9. (Color online) - Displayed quantities as well as the initial solution are the same as in Fig. 8. Now the initial noise is P = 0.7. 
Here the noise is strong enough to open a hole in the domain wall, which then emits ordinary vortices and decays by being absorbed by the 
boundary of the domain. 



